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Entrainment of randomly coupled oscillator networks by periodic external forcing applied to a 
subset of elements is numerically and analytically investigated. For a large class of interaction func- 
tions, we find that the entrainment window with a tongue shape becomes exponentially narrow for 
networks with higher hierarchical organization. However, the entrainment is significantly facilitated 
if the networks are directionally biased, i.e., closer to the feedforward networks. Furthermore, we 
show that the networks with high entrainment ability can be constructed by evolutionary opti- 
mization processes. The neural network structure of the master clock of the circadian rhythm in 
mammals is discussed from the viewpoint of our results. 
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INTRODUCTION 



The study of complex networks has applications in 
various fields including biology and engineering, and at- 
tracts growing attention 0, Q ■ In the last decade, much 
progress has been reached in understanding complexity 
of network architectures that represent skeletons of net- 
works. However, not only the architecture but also the 
dynamics taking place in a network is important 3] . This 
dynamics determines functions of the networks, related, 
e.g., to information processing in the brain (which is a 
network of neurons) or to the production process in a fac- 
tory (which is a network of machines). Most of the real- 
world networks are constructed in such a way that they 
implement certain desired dynamical behaviors. Not only 
the individual components of a network (such as single 
neurons), but also the network architecture (e.g., con- 
nections between neurons) should be appropriately de- 
signed. Thus, it should be expected that real-world net- 
works have architectures reflecting their desired dynam- 
ical behaviors. If one can extract topological properties 
of the networks relevant for their desired dynamics, this 
may provide insights into the meaning of network ar- 
chitecture from the viewpoint of dynamics and help to 
understand the design principles of functional networks. 

In the present study, the focus is on the synchroniza- 
tion of oscillators coupled into complex networks. Syn- 
chronizationplays a crucial role in functioning of various 
systems 0, IE IE S ■ O nc of the most intriguing ex- 
amples is the circadian (i.e., approximately daily) clock 
in mammals 0, ^| . The circadian rhythms of the whole 
body are orchestrated by a central clock in the brain, 
called the suprachiasmatic nucleus (SCN). This brain tis- 
sue is formed by a population of special neurons, known 
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as clock cells. A single clock cell exhibits a robust circa- 
dian rhythm in its firing rate and thus each such cell is a 
self-sustained oscillator (with oscillations determined by 
a regulatory loop in a group of genes at the single-cell 
level). In the SCN, neurons mutually synchronize their 
physiological rhythms even in absence of any environ- 
mental assistance ^lj, so that this cell population can 
generate collective periodic signal. Additionally, a differ- 
ent kind of synchronization takes place in the SCN which 
represents entrainment to environmental rhythms. (Here 
and below, we use the term "entrainment" to describe 
synchronization to an external input, distinguishing it 
from autonomous mutual synchronization). The entrain- 
ment of the SCN is essential for the proper functioning of 
this biological clock. Generally, the intrinsic period of the 
circadian rhythm is significantly different from 24 hours 
(in humans, it would typically be 25 h). Therefore, it 
must be tuned to the 24 h period through some external 
influences. Moreover, the phase of circadian oscillations 
needs to be locked appropriately to the local time. The 
entrainment to the environment is mediated by the light 
information coming from the eyes, which acts as external 
periodic forcing for the SCN. However, only a distinct 
subset of neurons receives and processes this light infor- 
mation, which implies that the rest of the population 
should be indirectly entrained via intercellular commu- 
nication (see Ref. j^] and references therein). There are 
various types of intercellular communication inside the 
SCN, including the communication through several types 
of neurotransmitters |lfj| | . The communication through 
neurotransmitters goes via chemical synapses of neurons 
forming a complex directed network. 

With respect to mutual synchronization of oscillators, 
the impact of the network architecture has been inten- 
sively investigated in recent years 0, 0, 0, 0, 0, 0, 
ITgj . However, only a few studies have dealt so far with 
the entrainment (i.e., synchronization to external peri- 
odic forcing) of oscillator networks 0, [53, • The ob- 
jective of the present paper is to identify principal topo- 
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logical properties that determine the entrainment ability 
of oscillator networks and to find the dependence of the 
entrainment ability on these topological properties. Al- 
though the problem is motivated by the SCN, our inter- 
est is more general. Therefore, we use the phase oscilla- 
tor model [f| , which provides a good approximation of a 
broad class of coupled oscillator systems and is expected 
to yield general insight into the entrainment phenomena. 
In the earlier Letter, a particular phase oscillator model 
has been studied and it was found that, in this model, 
the entrainment threshold for the coupling intensity be- 
tween oscillators increases exponentially with the depth 
of a network, defined as the mean forward distance from a 
pacemaker (i.e., source of external forcing) to the network 
nodes |2f)j |. Hence, it was found that the entrainment is 
practically only possible for the shallow networks with 
low hierarchical organization. Here, we give all details of 
the derivations, extend the analysis to different network 
models and demonstrate the universality of the earlier 
preliminary results. Additionally, a new class of random 
networks is considered, which allows further to explore 
the design principles of entrainable networks. The op- 
timization of the network architecture with respect to 
the entrainment through a dynamical learning process is 
further performed. 

The article is organized as follows: The model is intro- 
duced in Sec. ^ Then, we first show numerical results 
obtained for standard random networks in Sec. Ill II In 
Sec. IIVI the model is analytically investigated under a 
certain approximation for the network architecture. The 
comparison with numerical results is also provided. In 
Sec.0 a different class of random networks is introduced, 
i.e. directionally biased networks, and it is shown that the 
entrainment is strongly enhanced for the networks close 
to the feedforward- type. In Sec. I VII dynamical evolution 
of the network architecture using two kinds of learning 
algorithms is considered. The results are discussed, with 
a special emphasis on the design principles of biologi- 
cal clocks, in Sec. IVIII followed by main conclusions in 

sec. rvTm 

II. THE MODELS 

We consider a system of N + 1 phase oscillators, one 
of them being a pacemaker. The basic model is given by 
a set of evolution equations for the oscillator phases (f>i 
and the pacemaker phase </>o, 

JV 

= uj + — ]T AijT^i - (f>j) + fxBif^i - O ), 

<j) = lu + Auj. (1) 

The topology of network connections is determined by 
the adjacency matrix A whose elements Ay are cither 
1 or 0. The mean degree pN is the average number of 
incoming connections per node (and p is called the con- 
nectivity). The element with i = is special and rep- 



resents a pacemaker. Its frequency is increased by Auj 
with respect to the frequency u> of all other oscillators. 
The pacemaker is acting on the oscillators 1 < i < N\, 
the action being specified by the coefficients Bi taking 
1 for 1 < i < N\ and otherwise. The coupling be- 
tween elements inside the network is characterized by 
the coupling function T(x) and the (positive) coupling 
intensity coefficient n. In absence of a pacemaker, such 
networks undergo autonomous phase synchronization at 
the natural frequency ui if the coupling is attracting, i.e., 
if r'(0) < 0. The coupling to the pacemaker is character- 
ized by the 27r-periodic function T(x) and the (positive) 
intensity coefficient [i. We assume that functions Y(x) 
and T(x) are 27r-periodic, nonconstant, and smooth. 

Without loss of generality, our model can be simpli- 
fied. By going into a rotating frame, we have u> = 0. 
The maxima of the coupling function, maxT and maxT, 
are chosen equal to unity by properly defining the coef- 
ficients k and /i. Moreover, rescaled time t' = t Aui and 
rescaled coupling strengths k! = k/Aw,// = (i/Auj are 
introduced |30|. After that, the model takes the form of 
Eq. with Auj — 1 and to = (below, we drop primes 
in the notations for the rescaled quantities). Note that, 
in terms of the original model Q), an increase of the 
rescaled coupling between the elements corresponds ei- 
ther to an increase of coupling k or to a decrease of the 
relative pacemaker frequency Auj. 

The presence of a pacemaker imposes hierarchical or- 
ganization in the network architecture, which plays a cru- 
cial role in determining the entrainment ability. For any 
node i, its distance h with respect to the pacemaker is 
defined by the length of the minimum forward path sep- 
arating this node from the pacemaker. Thus, the whole 
network is divided into a set of shells, each of which is 
composed by oscillators with the distance h from the 
pacemaker. The shell population is given by the num- 
ber of the oscillators with distance h. The depth L of a 
network is introduced as 

L = ^Y. hN ^ ( 2 ) 

h 

which is the average distance from the pacemaker to the 
entire network. We may classify the types of connec- 
tions into forward, backward, and intra-shell connections, 
which are, respectively, connections from the nodes in a 
certain shell h to the nodes in the next shell h + 1, from 
the nodes in a certain shell h to the nodes in the shallower 
shells k < h, and between the nodes inside the same shell. 
We call the connections coming from a certain node the 
outgoing connections of the node, and the connections 
received by a certain node the incoming connections of 
the node. 

In the analysis of the model, several further assump- 
tions will be made. The number N\ of elements, directly 
connected to the pacemaker, is assumed to be small as 
compared with the total size N of the network. The mean 
degree pN is chosen large enough, so that the networks 
do not become disconnected. We assume T(0) — which 
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implies that coupling between connected oscillators van- 
ishes when their phases are synchronized. Moreover, cou- 
pling of the elements to the pacemaker is chosen to be 
much stronger than coupling between the elements in the 
rest of the network (fj. 3> k). In all the numerical simu- 
lations performed in the present paper, instead of using 
Eq. 0} wc put <pi = <fio for 1 < i < Ni, corresponding to 
the the limit /i — > oo. 



III. NUMERICAL INVESTIGATIONS FOR 
STANDARD RANDOM NETWORKS 

We begin our analysis of entrainmcnt phenomena by 
considering the important special case of standard ran- 
dom networks, also known as Erdos-Renyi (ER) networks 
[22I 123^ . These networks are generated by independently 
assigning with probability p for any pair i and j of 
the network nodes a connection that leads from the node 
i to the node j. Hence, elements of the adjacency 
matrix A are chosen to be 1 with probability p and 
otherwise, and A is asymmetric in general. Only sparse 
random networks with a small mean degree pN <C N will 
be considered. We use the following coupling function 



T(x) = 



sin(x + a) — sin a 



1 



sin a 



(3) 



where a is constant. Note that T(0) = and maxTfi) = 
1 for any a. Note also that for a = this coupling 
function is the same as that used in Ref. [20l |. 

Numerical simulations are started with random phases 
for the oscillators i > N\. For each oscillator, its effective 
long-time frequency u>i is computed as 



T 



(4) 



with sufficiently large T and to- 

Numerical simulations show that the response of a net- 
work to the introduction of a pacemaker depends on the 
strength k of coupling between the oscillators. When 
this coupling is sufficiently large, the pacemaker entrains 
the whole network (i.e., u>i = 1 for any i). Under en- 
trainment, the relative phases <pi — 4>q are locked. As the 
coupling strength n is decreased, the entrainment breaks 
down at a certain threshold value K cr (for a = 0, see 
Fig-0- Our simulations show that, in most cases, syn- 
chronization between the first and the second shells is 
the first to break down, and the frequencies of oscillators 
in the shells h > 2 remain equal for any k, i.e. u>i = ujj 
for i,j > N\. The results for a = 0.5 and a = —0.5 are 
similar. 

Figure [21 displays the thresholds K CT (in the logarith- 
mic scale) for a large set of networks with the fixed mean 
degree pN and different numbers N\ of oscillators in the 
first shell. As can be seen, the entrainmcnt threshold n CI 
may vary largely for different network realizations even 
if N\ is fixed. The results for a = —0.5, 0, and 0.5 are 
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FIG. 1: (color online). Long-time frequencies for various k. 
Data obtained for three samples of random networks is plotted 
with different symbols. The depths of the networks are L — 
3.01,3.33 and 3.45. The lines show the theoretical curves 
with k ci = 67.5, 125 and 169 respectively. TV = 100, p = 
0.1, Ni = 1, and a = 0. 



labeled as A, B and C. The same set of networks is used 
for all three values of a. For a given network, the en- 
trainment thresholds obtained for different values of a 
are close, implying that the entrainment threshold does 
not depend on the particular form of the coupling func- 
tion. Each group of networks with a certain Ni is dis- 
played with its own symbol. Each group is characterized 
by a distribution of depths L and generates therefore a 
cluster of data points. Correlation between the entrain- 
ment threshold K cr and the network depth L is evident 
even within a fixed Aq. The distributions inside each 
cluster and the accumulation of the clusters yield the 
dependence k C t:(L) of the entrainmcnt threshold on the 
network depth. Remarkably, the observed dependence is 
well numerically approximated by the exponential law 



K cr oc (1 +PN) 1 



(5) 



The same functional dependence is found for networks 
of other system sizes and different average mean degrees 
pN (see Fig. I3J) . The entrainment thresholds for N = 100 
and N — 200 can be well fitted by the same function 
with the same fitting parameter c, which suggests that 
the entrainment threshold depends not explicitly on the 
system size N, but on the average degree pN. 

Next we consider relaxation to the entrained state. 
The relaxation time for each generated network has been 
measured as follows. First, we have run a numerical sim- 
ulation for a long time starting with a random initial 
distribution of phases. If the relative phases tpi = <pi — cf>o 
for all i were not locked at the end of the simulation (im- 
plying that the network did not become entrained for the 
chosen coupling intensity), such a network has been dis- 
carded. Then, the simulation was repeated starting from 
the initial condition <f>i(t = 0) = ijjf + eyi for i > Ni, 
where yi is a random number independently taken from 
the uniform distribution within (0,1), e is a small coef- 
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FIG. 2: (color online). Dependences of the entrainment 
threshold on the depth L for a large set of random networks 
of size N = 100 and pN = 10. The coupling function is given 
by Eq. © with (A) a = 0.5, (B) a = 0, and (C) a = -0.5. To 
separate the data sets, here we have plotted vk ci with v — 5, 1 
and 0.2 respectively for (A), (B) and (C). The lines are the 
theoretical dependence cvk, ci , where k cv is given by Eq. (gSJ 
and c = 0.60 is a fitting parameter. 



ficient and ipf are the relative phases of oscillators in 
the entrained state. The time dependence of the dis- 
tance D{t) = J2i \/ {'tpiit) ~ V'i'} 2 nas been monitored. 
The relaxation time r was defined as the time t at which 
the ratio D(t)/D(Q) becomes equal to e _1 . Figure Q] dis- 
plays the relaxation times r (in the logarithmic scale), 
obtained numerically for a large set of networks and for 
a fixed coupling strength k. Again, exponential depen- 
dence on the depth is evident. This dependence is well 
fitted by the function 



oc (1+pN) 1 



(0) 



There is a divergence of the relaxation time around 
L = 3.7. This divergence occurs around the region of 
the depth where given coupling strength k is close to the 
entrainment threshold, i.e., n ~ K cr (L). Above this re- 
gion, entrainment rarely happens. 

Eigenvalues of an entrained state are displayed in 
Fig. [SJ They were obtained by numerically solving the 
stability matrix of our model QJ. We have assumed the 
limit \x — ► oo and perturbations are considered only in the 
subsystem i > Ni, so that there are N — N\ eigenvalues. 
The eigenvalue possessing the maximum real part is de- 
noted by A max . As seen, the magnitude of the real part of 
this eigenvalue, |ReA max |, is much smaller than those of 
others (which are distributed around \k\). We have pre- 
liminarily checked the eigenvectors and found that the 
associated eigenvector with A max corresponds to an ap- 
proximately identical phase shift of the whole subsystem 
i > Ni, and the rest eigenvectors correspond to relative 
motions inside this subsystem. This suggests the follow- 
ing relaxation process: Relative perturbations inside the 
subsystem quickly diminish within the characteristic time 
scale At -1 and, then, the phase differences between the os- 
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FIG. 3: (color online). Dependences of the entrainment 
threshold on the depth L for a large set of random networks. 
The coupling function is given by Eq. © where a — 0. (A) 
N = 200 and pN = 10 . The line is the the theoretical de- 
pendence ck ci . , where K cr is Eq. 133^ and c = 0.60 (the same 
as the value used for N = 100). (B) N = 100 and pN = 6, 20. 
The lines are the theoretical dependence CK cr , where 
Eq. 1331 and c = 0.67 and 0.45 respectively for pN = 6 and 
20. 



cillators in the subsystem become practically locked. The 
relaxation of the phase difference between the first shell 
and the subsystem proceeds slowly with the time scale 
| Re Amax | _1 - Such behavior has actually been observed 
in our numerical simulations. 

The results relating to the relaxation process do not 
change qualitatively for different values of a (we have 
further checked them for a = 0.5 and a = —0.5). Note 
also that the relaxation time does not depend on the 
particular form of small perturbations. 



IV. ANALYTICAL INVESTIGATIONS IN THE 
GLOBAL TREE APPROXIMATION 

In this Section, the model is analytically investigated 
using the heuristic global tree approximation. First, 
global tree networks will be introduced and then their 
similarity to the ER networks will be explained. After 
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FIG. 4: (color online). Dependence of the relaxation time 
on the depth L for an ensemble of random networks. The 
parameters are N = 100, p = 0.1, a — 0, and k = 300. The 
solid line is the theoretical dependence dr where r is Eq. 12811 
with K cr = 0.6pAT(l +pN) L ~ 2 (same as the theoretical lines in 
Fig- HJi an d d = 1.22 is an additional fitting parameter. The 
dotted line is the exponential function k ct /k oc (1 +pN) L . 
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FIG. 5: (color online). Eigenvalues (divided by n) of the 
entrained state. The parameters are N = 100, Ni = 2, p = 
0.1, a = 0, and k = 300 (the depth of the used network 
is L = 2.86 and the entrainment threshold for this network 
is found numerically to be K cr ~ 33). The real part of the 
maximum eigenvalue is Re(A max ) — 11. 



that, the analytical solution of the entrainment problem 
for global tree networks will be constructed. Note that 
this analytical solution will hold only in the limit when 
both pN and N are large. 



A. Global tree approximation 

A global tree network is a hierarchical network where, 
in any level, each oscillator has only one incoming con- 
nection from the higher level and exactly pN outgoing 
connections leading to the lower level. Additionally, each 
oscillator has exactly pN connections from the last (bot- 




FIG. 6: (color online). An example of a global tree net- 
work (JVi = l.pJV = 2,H = 4, and hence, N = 15). 
Thick solid lines, thin solid lines, and dotted lines repre- 
sents respectively forward connections, backward connections, 
and intra-shell connections. Numbers indicate the hierar- 
chical level of each node. This figure is drawn by Pajek 
jhttp: / /vlado.fmf.uni-lj .si /pub /netw orks /pajek/ ). 



torn) shell of this hierarchical network. No other connec- 
tions exist and the precise connection topology remains 
arbitrary. An example of a global tree network with 4 
hierarchical levels is shown in Fig. 

For such networks, shell populations grow exponen- 
tially, i.e. as Nh — Ni(pN) fl ~ 1 with the distance h from 
the network origin. The total number H of shells is de- 
termined by the condition Ylh=i ^ h = The depth 
of such a network is given by L = AT -1 Ylh=i = 
H + O(lfpN). Thus, for pN -> oo, the depth L coin- 
cides with the number of shells H. Therefore, H can be 
approximately replaced by L. 

Global tree networks share essential properties with 
the ER networks of large size TV ;§> 1 and high con- 
nectivity pN 3> 1. This similarity is briefly explained 
below (see Appendix ^] for a detailed discussion) . We 
first consider the pattern of forward connections of the 
ER network. By definition of our model, each node in 
the first shell receives one forward connection from the 
pacemaker. Each node typically gives pN outward con- 
nections. Thus, from the first shell, a total number of 
the outward connections is pNNi, on the average. All 
elements receiving a connection from the first shell form 
the second shell of the network. If the number of connec- 
tions leading from the first shell is much smaller than the 
total number N of elements to which they may lead (i.e., 
if pNNi <C N), each next outward connection from the 
first shell is typically received by a different node in the 
second shell. This means that, typically, an element in 
the second shell would be linked only to a single element 
in the first shell, as required by the tree structure. Con- 
sidering the third shell, we can notice that again that, if 
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FIG. 7: Schematic representation of the network structure 
along a forward path in the global tree approximation. Each 
circle with a number h is an oscillator of the shell h. Each 
right-arrow and each left-arrow represent respectively a for- 
ward connection and pN backward (or, intra-shell) connec- 
tions. 



the number N 3 = Ni(pN) 2 of out going; connections from 
this shell is small as compared with the network size TV, 
the tree structure would approximately hold for this shell 
too. 

Shell populations Nh grow exponentially with the num- 
ber h of the shell, i.e. N h = N^pN) 11 ' 1 . The tree struc- 
ture with respect to forward connections holds so far as 
these populations remain much smaller than the total 
network size. When pN is large, it can be shown that 
only the last two shells have populations of size O(N) 
and, thus, the tree structure approximately holds down 
to the third last shell. Analyzing further patterns of 
backward and intra-shell connections in a large ER net- 
work, we notice that populations of all shells, except the 
last of them, are of order o(N) and only the popula- 
tions of the two last shells are of order O(N). Therefore, 
each node outside of the last two shells receives back- 
ward connections mostly from these last two shells. On 
the average, the number of such connections is pN. Each 
node in the last two shells typically has pN incoming 
connections from the last two shells. This large number 
of connections between the last two shells strongly facil- 
itates synchronization of oscillators inside them, as can 
be seen in numerical simulations. Having this in mind, 
we may merge the last two shells in the network into a 
single shell. Once this additional (empirical) approxima- 
tion is used, the tree structure is extended to the entire 
network, including its last shells. 

Thus, in the global tree approximation the ER network 
is treated as a network having a. tree structure 1 with re- 
spect to forward connections and with the backward con- 
nections arriving only from the last network shell. In this 
network, every forward path starting from the pacemaker 
has a similar structure, as illustrated in Fig. [7| 

When this approximation is used, analytical solution 
can be easily constructed and its properties become more 
clear. Therefore, we prefer first to present our analy- 
sis using the global approximation, even also it remains 
partly empirical. The investigation for ER networks can 
also be performed without relying on the tree approxi- 
mation. While giving quantitatively better results, this 
other analytical investigation is however more compli- 
cated and, therefore, we give it separately in Appendix 
A. 



B. The entrained solution 

In the considered tree networks, our model (Q) can be 
written as 



ft = ^nrt <t> h k - 1 ) + ^ £ - (7) 

where h > 2; §\ is the phase of the oscillator i in the 
shell h; </>fc _1 is the oscillator giving forward connection 
to the oscillator i in the shell h; and the summation is 
taken over all the oscillator j in the last shell L. Be- 
cause each oscillator in any particular shell has the same 
pattern of connections, i.e., one forward connection from 
the previous shell and pN backward connections from the 
last shell. Thus, the entrained state of a network with 
phase synchronization inside every shell is possible. For 
such a state, the entrainmcnt solution problem is reduced 
into that in the one-dimensional oscillator array along a 
forward path starting from the pacemaker (see Fig. [7J. 
Putting fa = 1 (which is the entrainment condition) and 
<t>i = Oh for all the oscillator i in the shell h, we get the 
following algebraic equations: 

fit (01 - t) + kT(0i - e L ) = 1, (8) 



— T(6 h - 6 h -x) + KT(8 h -9 L ) = l for h > 2. (9) 
pl\ 

For large pN, we may linearize T{0h—0k) as T'{<S){6h— Ok) 
for h, k > 2 [it will be shown that 82 — 0l is at most of 
order 0(l/pN) in the solution under entrainment, so that 
this linearizion is justified]. We may then solve Eq. @ 
from the last shell upward. For h = L, we get 

r(0r - e L -x) = K-^N. (10) 

Substitution this into Eq. © for h = L — 1 results in 

r(0 £ _i -0L-2) = K^pNii+pN). (11) 

From the following shells, we need to employ the linicar- 
lization, i.e., we approximate 

L 

r(9 h -e L ) = - £ T(6 k+1 -e k ) ^h>2. (12) 

k=h-l 

We then obtain 

r(6 h - 0fc_i) = «T VV(1 + pN) L - h for h>2. (13) 

Substituting Eq. (JT3J for h = 2 into Eq. ©, we also 
obtain 

(JL H 
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Note that the explicit expression for T(0i — t) is not 
needed. 

The existence conditions of the entrained solution are 
f (01 - t) < 1 and r(6 h - 9 h -i) < 1 for h > 2. The 
former condition is satisfied for p, 3> K (which is our 
assumption). Among the terms T(6h — @h-i), the term 
T(02 — 0i) is the largest one. The solution thus exists if 



r(0 2 - 



< 1. 



(15) 



Therefore, the entrained state of the network is possible 
only if the coupling intensity k satisfies condition k > K cr , 
where 



K cr = P N(l+pN) 



L-2 



(16) 



The result 1)16(1 is remarkable. According to it, the 
coupling threshold 

IS determined exclusively by topo- 
logical properties and does not depend on a particular 
form of the coupling function. Moreover, it is given only 
through a combination of pN and L. This fact suggests 
that the essential parameters of the system are pN and L 
rather than N, p and N± , in agreement with the previous 
numerical results obtained for the ER networks. 

From Eq. (fTfi|t . K C r for our tree network under consid- 
eration is roughly estimated as 



k c , ~ (pN) 



L-l 



N/Ni 



(17) 



Thus, large number for small Ni (as we assumed). 

This property will be used in the derivation of the relax- 
ation time. 

One can verify that 02 — &l is at most of order (pN)^ 1 
in the entrained state by substituting K cr into k in 
Eq. H13|) . Note also that although the phase differences 
&h — 0/i+i arc small in deeper shells, they should not be 
neglected when the cntrainmcnt threshold is derived, be- 
cause the phase difference 0i — 02 (which determines the 
existence condition of the entrained solution) is a con- 
sequence of the exponential growth of such small phase 
differences from the deepest shell upwards. 



C. Stability analysis of the entrained state 

To analyze stability of the entrained solution, only the 
limit p, — ► 00 is considered (so that 0, = 0i — > t for 
i < Ni). In numerical simulations, we have observed that 
the relaxation process is characterized by two distinct 
time scales: the first one characterizes fast relaxation in- 
side the subsystem i > N\ and the other corresponds to 
slow relaxation in the phase difference between the first 
shell and the subsystem. In this section, we first investi- 
gate the internal stability of the subsystem. Then, taking 
advantage of the separation of time scales, we hcuristi- 
cally construct the effective dynamical equation of the 
whole subsystem. 

We consider small perturbations for phases in the en- 
trained state, i.e., <$? = &h + $4% where Scj)^ is small 



perturbation. Because of the assumption p, — > 00, we 
put 5(j>\ = 0, i.e., $ = 0i = t. The model Q for h > 2 
then gives 

Hi = ^ 7 r'(0,-0,_ 1 )(^-^ 1 ) 



Sift). (18) 



Because 0/, — 0l = 0(l/pN), we approximate T'(9h — 
6 L ) ~ r'(0). In addition, because of the large number pN 
of connections from the last shell, we may approximate 



N L 

j 3 
Using these two approximations, Eq. I jlSjl reduces to 



(19) 



pN 



T'(6 h -0 h _i) {6$ - 8<t> h k - 1 ) + kV (0) {5$ -5ct> L ). 

(20) 



Since the first term is negligibly small for large pN for 
any shell h > 2, Eq. (|20|l is further approximated as 



(21) 



Thus, in our approximation, all eigenvalues associated 
with the relative motion inside the subsystem are degen- 
erated into 



A = «r'(o), 



(22) 



with degeneracy N— Ni — 1 . The dynamical equation 1(21(1 
describes the fast relaxation inside the subsystem: oscil- 
lators quickly relaxes to the average perturbation 6<f> L in 



the last shell, i.e., 



b L . After this fast re- 



laxation, the phase differences inside the subsystem are 
almost locked. 

Now we heuristically construct the dynamical equation 
for the phase difference between the first shell and the 
subsystem. We take the coordinate of the subsystem on 
its surface, i.e., the phase of the second shell, defined as 
^ = 02 + S(f) L . Because the total external force applied to 
the subsystem is given by N 2 kT(^ - t)/pN = NikT(^ - 
t), the effective dynamical equation for the surface ^ 
reads 



iVe* = JVi«;r(*-t). 



(23) 



Here, the effective size N e of the whole network is deter- 
mined by the condition that ^ = 02 when the subsystem 
is entrained (i.e. ^> — 1). From this condition, comparing 
Eq. (O for h = 2 and Eq. J23JI, we obtain 



N c = N lKc 



(24) 



Thus, from Eqs. (f^ and (|2"I|) the last eigenvalue is found 
to be 



A„ 



KT'(02-0l) 



(25) 
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The eigenvector of A max corresponds to the identical 
phase shift for all the oscillators i > Ni. 

For the phase difference 62 — 0i, there are at least 
one pair of steady solutions, one of which always sat- 
isfies r'(0 2 - 61) < 0. Thus, provided that r'(0) < 0, 
a stable entrained solution exists for k < k ct . This im- 
plies that the entrainment breakdown occurs only via the 
disappearance of the solution at n = K cr . 

Because large number for small N\ [see 

Eq. ltT7|) ]. the time scales are well separated, i.e., 
|A max | <C |A|. This fact justifies the approximation em- 
ployed in this subsection. In addition, the relaxation time 
against any general perturbation is thus given simply by 



(26) 



kT'{6 2 - 0i) 



Hence, in general, the relaxation time depends on the ex- 
plicit form of the coupling function, being different from 
the case of the entrainment threshold. Nevertheless, the 
dependence on the L is approximately the same as k ct 
in the region of L where r'(02 — 61) does not vary much 
with L. 

For the special case T(x) = — sin a;, Eq. H'25J) reduces 
into 



V k2 K cr 



The relaxation time is thus 



V^ 2 



(27) 



(28) 



Because K cr increases with L, n cr eventually coincides 
with k at certain critical L, which results in the diver- 
gence of r. In other region of L, the dependence of the 
relaxation time on the depth L is approximately the same 
as K cr , i.e. exponential. 



D. Below the entrainment threshold 

Dynamical behavior just below the entrainment 
threshold is considered. We again choose the limit fi — + 
00. Because of the property |A max | <C |A|, we only need 
to consider the dynamics of the subsystem, Eq. I|23|) . We 
introduce the slow mode x = VP — t and the bifurcation 
parameter e = (n — K cr )/ K cr- Substituting them into 
Eq. we obtain 



x = (1 +e)T(x) - 1. 



(29) 



Expansion of T(x) up to the second order around its max- 
imum x = x max yields 



r". 



■(x — x max ) 2 + higher orders, 



(30) 



where rjj' lax < is the second derivative of T(x) at 
x = x max . The steady solution of x (corresponding to 



the entrainment) thus disappears via a saddle-node bi- 
furcation at e = 0, and x oscillates with the negative 
frequency w x for e < 0. This frequency is found through 
the integration of Eq. (f 301) : 



2tt 



dx 



V" 

max 



(x 



c)72 



+ 0(e), (31) 



followed by lo x = — \/2e/r(^ iax . Hence, the effective fre- 
quency LOi of each oscillator in the subsystem (i > Ni) 
for e < is 



f 2e 
T" 

max 



0(e). 



(32) 



E. Comparison with numerical results 

The analytical results obtained under the global tree 
approximation are compared with the numerical results 
obtained for the ER networks. The theoretical curves 
plotted in Fig. ^ correspond to the function (|32|) for ap- 
propriate values of K C r, showing excellent agreement with 
the numerical data for n close to « cr , as expected. In 
Figs. |2 and |3] we display 



cpN(l +pN) 



L~2 



(33) 



where c is a fitting parameter added to our theoretical de- 
pendence qiufl. As already noticed in Ref. [2(j, although 
the principal dependence on the depth, (1 + pN) L , is 
correctly reproduced by the global tree approximation, 
the coefficient is somewhat different. We have thus in- 
troduced the fitting parameter c. As seen in Figs. and 
13 the principal dependence on the depth obtained for 
the tree network agrees excellently with the numerical 
results obtained for the ER networks. Note that the the- 
oretical dependence, analytically derived directly for the 
ER network in Appendix^ agrees well with numerical 
data without any fitting parameter. For the relaxation 
time, taking into account the correction for /t cr , we use 
instead of Eq. (|28|l the following function 



(34) 



y^K 2 ~ (cK CJ: y 



where the value of c is obtained by numerical fitting 
for the entrainment threshold and d is an additional fit- 
ting parameter. As seen in Fig. 01 the function (|34J) fits 
very well to numerical data. In Fig. El real parts of the 
eigenvalues except A max are scattered around k, as ex- 
pected from Eq. (J23) [note that T'(0) = 1 for a = 0]. 
Putting k = 300 and K cr = 33 into Eq. l(57|) . we obtain 



A„ 



-9.0. which is close to the numerical result shown 



in Fig.E|(A r 



-11). 
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V. ENTRAINMENT IN DIRECTIONALITY 
BIASED NETWORKS 

We can understand by looking at Eq. why entrain- 
ment is difficult in the networks with larger depths. The 
first term on the left side in the equation describes the 
force of the forward connection. Its sign is positive, so 
that it contributes to the entrainment. On the other 
hand, the second term describes the force from the back- 
ward connections, whose sign is negative. By moving 
this second term to the right side, it is seen that it es- 
sentially increases the frequency to which the oscillators 
synchronize. In other words, the backward connections 
act as a load for the entrainment. Moreover, the number 
(pN) of the backward connections is much larger than 
the number (one) of the forward connections. Thus, to 
compensate such strong unbalance, the phase difference 
associated with the forward connection {Qh-i — Oh) needs 
to be much larger that that associated with the backward 
connections (0^ — #£,). The effect accumulates exponen- 
tially along a forward path of the length L starting from 
the pacemaker and ending at the last shell. This ac- 
cumulation results in the exponential growth the phase 
difference from the last shell upwards, and thus, the ex- 
ponential dependence of the entrainment threshold on 
the depth L. Hence, due to the large number of back- 
ward connections, the entrainment is very difficult for the 
networks with large depths. 

It is thus expected that the entrainment threshold de- 
creases significantly for the networks that are closer to 
the feedforward architecture (i.e., to a network without 
backward connections) . This is indeed demonstrated be- 
low using a special class of random networks which we 
call directionally biased networks. To construct a direc- 
tionally biased network, we first generate a directed ER 
random network and choose N\ nodes as the first shell. 
Then, we redefine its backward connections. Namely, for 
every backward connection, we decide to retain it with 
probability £ or delete otherwise. We call £ the backward 
connectivity of the network. Note that a directed ER net- 
work corresponds to £ = 1 and the feedforward network 
is obtained for £ = 0. 

To solve the entrainment problems for such networks, 
the same approximations as in Sec. HVI are applied. The 
global tree approximation for the forward connection pat- 
tern can be used. Furthermore, it can be assumed that 
all backward connections come from the last shell and 
that the number of such connections received by a os- 
cillator is £pN. The former approximation is applicable 
when pN is large and the latter is applicable when £pN 
is large. Thus, we require that £pN is large. Under this 
condition, the same linear approximation as in Sec. IIVI 
can be used, because the phase difference 02 — 0l turns 
out to be at most of order (^pN)^ 1 . The entrained so- 
lution for the approximated network is found by solving 
algebraic equations 



— T(0 h - 9h-i) + i^T(9 h -9 L ) = l for h > 2. (36) 
pN 

We obtain 

Ftfh-Oh-i) = K- 1 pN(l + &N) L - h hih>2, (37) 



= pN(l + £, P N) 



L-2 



(38) 



The stability analysis is performed in the same manner, 
leading to 



(39) 



kT' 



/if (0! - 1) + <r(0! - 9 L ) = i, 



(35) 



Hence, both the entrainment threshold K cr and the re- 
laxation time r decrease dramatically as the backward 
connectivity £ gets smaller. The role played by £ is more 
significant for the networks with a higher hierarchical or- 
ganization (i.e., with a larger depth L). 

The employed linear approximation is applicable only 
for large £pN . For £ = 0, the model is however ex- 
actly solvable without the linear approximation. Since 
there are no backward connections in this case, perfect 
phase synchronization inside each shell is possible. Con- 
sequently, we get Eqs. (TT7JI and also for £ = 0. Triv- 
ially, the phase difference Oh-i — Oh = A between the 
neighboring shell is constant for h > 2, given by the re- 
lation r(— A) = pN/n. The stability analysis of this so- 
lution is straightforward. The eigenvalues are ^,T'(0i —t) 
(multiplicity N x ) and KT'(-A)/pN (multiplicity N-Ni). 
For k > K cr , there is at least one pair of solutions, one of 
which satisfies r'(— A) < 0. Therefore, for any T, there is 
always a stable solution for k > k ci . Note that this state 
has a constant phase gradient and can thus be described 
traveling wave. 

Numerical results for directionally biased random net- 
works arc displayed in Fig.|S| For the entrainment thresh- 
olds (Fig. the theoretical dependence (1 + £,pN) L 
fits well the numerical data. For the relaxation time 
(Fig. |HJ3), the theoretical dependence (1 + £,pN) L also 
fits well the numerical data, except for £ = where a lin- 
ear function fits best. This linear dependence is natural 
because the solution is wave-like and it is thus expected 
that the time necessarily to transmit information is pro- 
portional to the system length (i.e., to the depth L). Note 
that the entrainment breakdown occurs near L = 4 for 

Similar dependences can be obtained for the weighted 
networks. To construct them, we first generate an ER 
network and introduce connections from the pacemaker. 
We then put Aij = £ for all existing backward con- 
nections, while keeping the weights of other connections 
equal to cither or 1. For such networks, we also ob- 
tain Eqs. I|37|) - (|39|) because the same approximations are 
applicable. 
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FIG. 8: (color online). Dependences of (A) the entrainment 
threshold and (B) the relaxation time on the depth for an 
ensemble of directionally biased networks. The solid lines are 
the exponential functions proportional to (1 + ^pN) L . The 
dash line is a linear function aL + b with appropriate fitting 
parameters a and b. Other parameters are k = 600, pN = 10, 
(A)N = 100 and (B)N = 400. . 



VI. EVOLUTIONARY OPTIMIZATION 

So far, the relationship between the entrainment 
threshold and the network depth has been established 
only for the special kinds of random networks. There- 
fore, it is natural to ask whether it also holds for the 
arbitrary complex networks. Obviously, we cannot an- 
swer it by investigate all possible types of complex net- 
works because their number is too large. Instead, an 
alternative way is chosen: we shall construct networks 
with lower entrainment thresholds (and thus better en- 
trainment ability) through an optimization process start- 
ing from an arbitrary network architecture and analyze 
changes of their topological properties in the course of 
evolution. If these networks actually evolve towards be- 
ing shallower and closer to a feedforward type, one can 
conclude that both properties are generally essential for 
good entrainment ability. 

During the evolution, several topological properties 
shall be monitored. One of them is the depth L, which 
has already been defined. The other property is the back- 



ward connectivity £ defined as the ratio of the total num- 
ber of backward connections at each trial and that of the 
initial network. For comparison, we also introduce the 
forward connectivity Xi given by the ratio between the 
total number of forward connections in each trial and 
that of the initial network. The total number of outgo- 
ing connections from the first shell, denoted by n out , shall 
also be monitored. 

Two types of optimization algorithms will be cm- 
ployed. The first one is the straightforward optimization 
where the network structure evolves in such a way that 
the mean frequency of the whole network becomes closer 
to that of the pacemaker. In the second algorithm, each 
oscillator selects its incoming connections in such a way 
that its own frequency becomes higher. It will be shown 
that, for both kinds, the evolving network improves its 
ability to become entrained and, at the same time, indeed 
becomes less hierarchical and closer to feedforward net- 
works. Note that, concerning biological clocks, the first 
and the second algorithms could imitate respectively the 
evolution process of the neural network of the SCN for 
a species and the growing process of the neural network 
for an individual. A more biological algorithm (or, learn- 
ing process) has been employed elsewhere and resulted 
in similar results |24j . 

The initial setup is common for both types. The pro- 
cess starts from a random ER network. We put fa (t) = t 
for i < Ni. Throughout the evolution process, the con- 
nection pattern from the pacemaker is maintained, so 
that the first shell is always composed of the oscillators 
i < Ni. The coupling strength k is fixed far below the 
entrainment threshold for the initial network. Thus, all 
the oscillators except in the first shell initially have fre- 
quencies close to their natural ones. 

The first evolution algorithm. A numerical simulation 
is run starting from random initial phases. The long- 
time frequency (lu) = (TN)' 1 £\ [fa(2T) - fa(T)] , aver- 
aged over the whole population and for sufficiently long 
time T, is determined. Then, the adjacency matrix A 
is mutated. This is done by "rewiring" : Wc choose ran- 
domly an existing directed link (Aij = 1) and eliminate 
it (A^ — > 0), and then, choose randomly a missing link 
(Ai'ji = where i' ^ j') and add a new connection there 
(Aiiji — > 1). After the mutation, a numerical simulation 
is again run starting from new random initial phases and 
the new frequency (lu)' is determined. If (lu)' is closer 
to 1 than (lu), the mutation is accepted. Otherwise, the 
network is reset to its structure before the mutation. The 
iteration process is repeated until the long-time frequency 
(lu) becomes equal to unity. 

A typical dependence of the average frequency and 
two different topological properties during an evolution 
is shown in Fig. [3J\.. The tendency to decrease both the 
depth L and the backward connectivity £ is evident. In 
the most of evolutions where the average frequency (lu) 
has increased, a decrease in L and/or £ was observed. 
In Fig. |5j3, evolution of two other topological properties 
is shown. We see that the correlation between the feed- 
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ing incoming link of this oscillator (A^ = 1), delete it 
(Aij — > 0), then choose randomly a missing link to this 
oscillator (Aiji = where j' ^ i) and add a new link 
there (Aiji —* 1). After the mutation, a numerical simu- 
lation is repeated and the new frequency of the oscil- 
lator i is measured. The mutation is accepted if w[ > tt>,, 
and rejected otherwise. At the next step, we again ran- 
domly choose an oscillator and repeat the same proce- 
dure. Thus, in this evolutionary process, the mutation 
is done according to the individual activities of the os- 
cillators. Note also that, in this evolution algorithm, the 
total number of incoming connections for each oscillator 
is maintained constant. 
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FIG. 9: (color online). Evolution process under the first al- 
gorithm. Parameter values are N = 100, p = 0.1, Ni — 1 and 
k = 20. The depth of the initial random network is L = 3.43. 
(A)The solid line is the long-time frequency (ui) averaged over 
the whole network. The dash line and the dot line are respec- 
tively the backward connectivity £ and the depth L. (B)Two 
other topological quantities are plotted. The solid line and 
the dot line are respectively the total number n ou t of outgoing 
connections from the first shell and the forward connectivity 
X- 



forward connectivity x an d the mean frequency (u>) is 
weak. An increase in n ou t indicates the emergence of a 
hub (i.e., of a node with a large number of outgoing con- 
nections), which strongly contributes to decreasing the 
depth L. It is worth noticing that an increase in (w) does 
not necessarily imply an increase in n on t- Apparently, the 
entrainment ability relies on more fine properties in the 
network architecture, best characterized through L and 
£. In our numerical simulations, we have tried several 
sets of parameter values and several different initial ran- 
dom networks, always obtaining qualitatively the same 
results. 

The second evolution algorithm. For each iteration 
step, an oscillator i is randomly chosen. A numerical 
simulation is run and the long-time frequency of this os- 
cillator, Wi = [<f>i( 2T ) " <M T )]/ T with sufficiently large 
T, is determined. Then, a structural mutation of the 
network is introduced. We choose randomly one exist- 



Figure ITUR displays typical dependence of the aver- 
age frequency and two topological properties of the net- 
works during a single evolution. Although the average 
frequency (uj) of the whole population does not always 
increase in each iteration step, the network architecture 
changes during the evolution similar to what has been 
found for the first algorithm, towards the networks of 
smaller depth L and smaller backward connectivity £. In 
Fig. HUb. emergence of a hub and weak correlation be- 
tween (u>) and x are seen, similarly to the results in the 
first algorithm. It should be emphasized that, in the sec- 
ond algorithm, globally coordinated network architecture 
emerges solely through the "judgments" of individual os- 
cillators based only on their own activity. The reason 
why this local optimization rule leads to the globally co- 
ordinated structure is that the oscillators are mutually 
frequency-synchronized in most cases and therefore their 
individual frequencies usually coincide with the average 
frequency of the whole population. Hence, the second 
algorithm also works similarly to the first one. 



In both types of algorithms, our numerical results show 
that the development of the entrainment ability is fol- 
lowed by decreases in the depth L and the backward con- 
nectivity £, suggesting that these two topological prop- 
erties are the primary factors in determining the entrain- 
ment ability. It can be also noticed that the evolving 
networks tend to develop strong heterogeneity in their 
outgoing degrees. The depth of an evolving network be- 
comes smaller via an increase in the outgoing degrees of 
nodes in the shallow shells. In this manner, hubs of out- 
going connections are formed. On the other hand, the 
development of small backward connectivity is mainly 
due to a decrease of the outgoing degrees in the deep 
shells. Consequently, the heterogeneity of outgoing de- 
grees tends to be larger as the optimization process goes 
on|3lj. Many real- world networks, including scale-free 
networks, are known to have strong heterogeneity in their 
degrees Q. Our results suggest that evolving networks 
with good entrainment ability also become strongly het- 
erogeneous. 
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FIG. 10: (color online). Evolution process under the second 
algorithm. The depth of the initial random network is L = 
3.22, corresponding to k ci ~ 100. Parameter values and the 
definitions of the exhibited lines are the same as in Fig. |U] 



VII. DISCUSSION OF BIOLOGICAL ASPECTS 

As already mentioned in Sec. [I] our study is moti- 
vated by the SCN, the tissue orchestrating the circadian 
rhythms of the whole body in mammals. In this section, 
we discuss the neural network structure of the SCN and 
its possible roles from the viewpoint of the results of the 
present investigations. 

The SCN is anatomically organized into two groups, 
often called the "core" subdivision and the "shell" sub- 
division (see, e.g., Ref. ^3)- Only the "core" subdivision 
receives photic input and, thus, this part corresponds to 
the first shell in our model. The "shell" subdivision corre- 
sponds to the rest of the shells in our model. Within and 
between the "core" and the "shell" subdivisions, there 
are several types of inter-cellular communications which 
influence circadian oscillations |1C( . The pathways (i.e., 
the network structure) of each type of communication 
may be different. 

From the viewpoint of our results, if the SCN is con- 
structed so as to reach the best ability for entrain- 
ment, the network architecture should be feedforward. 
This means that the unidirectional connectivity from the 



"core" subdivision to the "shell" subdivision should be 
found. This is indeed the case for communication via 
VIP (vasoactive intestinal polypeptide), which is a neuro- 
transmitter released only by neurons in the core subdivi- 
sion. VIP is one of the leading candidate factors for the 
synchronization pathway inside the SCN |25j . However, 
VIP is not the only one communication agent. Com- 
munication via GABA (7 -aminobutyric acid) is further 
possible for almost all neurons in the SCN and seems to 
play a crucial role in achieving synchronization between 
the "core" and the "shell" subdivisions [H. The GABA 
connection pattern is not feedforward but bidirectional 
between these subdivisions. However, remarkably, it has 
been conjectured from experimental studies that coupling 
from the "shell" to the "core" is weaker than in the other 
direction [2(|, suggesting that the network would effec- 
tively be close to feedforward one also with respect to the 
GABA communication. 

It is known moreover that, although the response of the 
"core" structure to a sudden phase shift in the environ- 
ment is very fast, the response of the "shell" subdivision 
to such phase shifts is significantly slower [26t l27j . From 
the viewpoint of our results, it can be conjectured that 
the slow response of the " shell" subdivision is due to the 
more hierarchical, non-feedforward organization of neu- 
rons inside this subdivision. As we have seen, if a network 
is hierarchical, the response can become very slow even if 
the first shell is strongly connected to the environment. 

It is interesting to ask why the "shell" subdivision, 
coupled to the "core" , actually exists in the SCN and why 
this subdivision could be more hierarchically organized. 
One reason may be that such an organization is needed to 
keep the autonomy of the biological clock. The SCN must 
be capable to synchronize to the environmental rhythm, 
but, on the other hand, should not be too sensitive to 
the environmental information. In other words, the SCN 
must possess both autonomy and adaptivity in a good 
balance. The hierarchical organization and directionality 
of the network architecture in the SCN could be useful 
to reconcile these contradicting requirements. 



VIII. CONCLUSIONS 

Our main result for the random ER networks and the 
random directionally biased networks is schematically il- 
lustrated in Fig. 1111 where k is non-rescaled coupling 
strength, i.e., k in the original model Q). The cntrain- 
ment occurs in a gray region. For the ER networks, this 
region becomes exponentially smaller for the networks of 
higher hierarchical organizations (i.e., with larger depth 
L). Thus, in practice, the entrainment is possible only 
for shallow networks. For the directionally biased net- 
works, the entrainment region becomes significantly en- 
larged for the smaller backward connectivity £ even if the 
networks are hierarchical. The feedforward network is 
the best one for the entrainment. The relaxation time to 
the entrainment has approximately the same dependence 
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FIG. 11: Schematic representation of the entrainment win- 
dow, inside which the entrainment occurs. The slopes of the 
edges of the entrainment window for u > and uj < are 
proportional to maxF and minF respectively. 



on the network architecture as that of the entrainment 
threshold. These results are general and hold for a large 
class of coupled phase-oscillator systems (and thus also 
of weakly coupled limit-cycle oscillator systems) with at- 
tracting couplings. The networks with such topological 
properties are shown to emerge naturally through differ- 
ent kinds of evolutions aimed at increasing the entrain- 
ment ability. 
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APPENDIX A: ANALYTIC INVESTIGATIONS 
WITHOUT THE GLOBAL TREE 
APPROXIMATION 

In Sec. IIVI to solve the entrainment problems we have 
used the heuristic global tree approximation for the ran- 
dom ER networks. By a different method, entrainment 
problems can be solved directly for the ER network of 
large size N 3> 1 and large degree pN ^> 1. Although 
this systematic derivation is more accurate, it is also more 
technical and, therefore, we have preferred to present 
first, in the main part of the article, the analysis based 
on the global tree approximation. The results of this 
derivation agree with what we have found before in the 
global tree approximation and can be viewed as providing 
further support for it. 

We begin by estimating the number muk of incoming 
connections leading from all nodes in the fcth shell to a 
node in the hth shell. The network pattern of forward 
connections is first considered, i.e., rrih^h-i for ft. > 1 is 
estimated. By definition, there are Ni nodes in the first 
shell, each of which receives only one forward connection 
(coming from the pacemaker). We thus have mifl = 1- 



The expected total number of outgoing connections from 
the nodes in the first shell is NipN because each node 
gives typically pN outgoing connections. If pNNi <C N, 
every outgoing connection almost surely connects to each 
individually different node outside the first shell. There- 
fore, in a good approximation, the number N 2 of nodes 
in the second shell is NipN, and each node in the second 
shell receives only one forward connection, i.e., r/12,1 = 1. 
The same property holds up to a certain shell ft* , where 
for the first time the total number of outgoing connec- 
tions from the shell becomes of the order O(N) or larger, 
i.e., pNN h *-i = o(N) and pNN h , > 0(N). We thus 
have N h = N 1 (pN) h - 1 (<^ N) and m h ,h-i = 1 for 
ft, < ft* . In other words, the network pattern of forward 
connections takes a tree structure in a good approxima- 
tion up to the shell ft*. The population Nh*+i of the 
next shell is of the order of O(N). Since J2h<h* is of 
o(N), almost all outgoing connections from the shell ft* 
connect to nodes in the shell ft,* + 1. Therefore, the ex- 
pected number of incoming connections from the shell ft* 
to a node in the shell ft* + l is pNNh* /Nh*+i = fhh*+i,h* ■ 
The statistical deviation from this expected number may 
not be neglected because fhh*+i,h* is not a large number 
in general. Since Nh*+\ = O(N), every node in the net- 
work almost surely receives connections coming from the 
shell ft* + 1, which implies that the rest of nodes belong 
to the shell ft* + 2. The expected number fhh*+2,h*+i of 
forward connections to a node in the last shell is pNh* + \. 
This is a large number of 0(pN), and thus the statis- 
tical derivation from the number can be neglected, i.e., 
m h * + 2.h' + i =pNh*+i- 

Since L = E h hN h/N ~ {(ft* + l)N h . +1 + (ft* + 
2)N h * +2 }/N ~ ft* + 1 + N h , +2 /N, we get 



[L] 1, 



(Al) 



where [L] is the integer part of L. According to the result 
in [2^ , the depth of the ER random network with large 
size N is estimated as 



L 



]n(N/Ni) -7 
\n(pN) 



1.5, 



(A2) 



where 7 ~ 0.5772 is the Euler constant. 

Next, the structure of backward and intra-shcll con- 
nections is considered. A node in the shell h typically re- 
ceives pNk connections from nodes in the shell k (k > h). 
The numbers of such connections from the last two shells 
are of 0(pN). Relative statistical deviations from these 
numbers are of order (pN)^ 1 ^ 2 and thus negligible. We 
thus obtain m/^. = pNk for fc = h* + 1, h* + 2 and h < k. 
The number of backward connections from other shells 
is of o(pN), and thus negligible, i.e., we approximate 
mh,k = for h < k < h* . 

Our estimation for the network structure is sum- 
marized as follows. The shell populations Nh are 
N^pN) 11 - 1 for 1 < h < h* and of order 0(N) for 
h = h* + 1 and h* + 2. The number m^. & of incom- 
ing connections of a node in the shell h from the nodes 
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in the shell k is given by rrihji-i — 1 for 1 < h < h* and 
by m h , +1Jl , = pNN h ,/N h , +1 ; m h . +2 ,h*+i = pN h * +1 ; 
rn-hk = for k < h — 1 (by definition); mh,k = pN^ for 
k = h*+l, h*+2 and h < k; and m h ^ = for h < k < h*. 

Now we solve the entrainment problem for the net- 
work under consideration. We assume that the phase 



difference 



between any pair of oscillators in 



the shells h > 2 is so small that we may linearize as 
T((j>i — 4>j) = r'(0)(</>i — <j)j). This assumption is justified 
later [it will be shown that \ fa — <f>j\ < 0(l/pN)]. 

Dynamics in the shells h > 1 is considered. We denote 
the average phase of the oscillators inside the last two 
shells h* + 1 and h* + 2 by 9 e d ge - Because the number 
of connections from the last two shells is large, we may 
approximate the coupling from those shells as 



1 

pN 



J2 AijT'iO)^ - fa) ~ r'(0)(^ - odgc ), (A3) 



where j denotes the oscillators inside the last two shells. 
Since every oscillator in the shell h < h* has effectively 
the same number of incoming connections from each in- 
dividual shell, a state with phase synchronization inside 
each shell h < h* exists. We denote the phase of the 
oscillators in the shell h by Oh- Under entrainment (i.e., 
4>i = 1), the phases of such a state can be found 
solution of MlgcbrMic equMtions 

JLV(O)(0 h - 6 h ^) + KT'(0)(9 h - cdgc ) = 1, (A4) 

where 2 < h < h* . The phase </>° ds ° of the oscillator i 
inside the last two shells h* + 1 and h* + 2 are found by 

-^^^r'(o)(^ dsc -^f)+ K r'(o)(0f so -0 cdgo ) = i. 

j 

(A5) 

Averaging Eq. (|A5() for all the oscillators inside the last 
two shells, we get 



^r'(o)(0 edge -M = i- 



(A6) 



Note that because of the strong internal coupling inside 
the last two shells, the oscillators inside the last two shells 
are almost phase-synchronized, i.e, <p^ ~ # e dge for h = 
h* + 1 and h* + 2. From Eq. (jA4j and Eq. {SB, for 
2 < h < h* we obtain 

N 



cdgc h ' ~ T'(p)KN h .' 

{1+pN) 



P N t-x , ^h*-h 



r'(o)«; 

pN 



N 

~Nh* 



T(02-0i) = —{1+pN) 

K 



h*-2 



1 + 



The existence condition of the entrained solution arc 
f(0i - 0o) < 1 and r(6» 2 - 6 X ) < 1. The former is sat- 
isfied for sufficiently large fi (which is our assumption). 



Thus, for a given network, the existence condition for the 
entrained solution is k < K rr , whore 



k ci = a - l pN(l+pN) h '- 2 (1 + ^t- 



N 
A 7 ^ 



(A8) 



Substituting K cr into k in Eqs. one can confirm 

that the phase difference 82 — Oedge is actually at most 
of 0(l/pN), so that the linear approximation is justified 
for large pN. 

In Eq. HA8jl . because h* increases with the depth L, 
the entrainment threshold approximately increases ex- 
ponentially with the depth L, although the effect of the 
last term l + N/Nh* is unclear. We thus employ further 
approximations and express the entrainment threshold 
(IA8|) as a function of the depth L and the mean de- 
gree pN. From Eq. (|A2|) . the network size can be ex- 
pressed as N — e 7 A r i^/V) L-1 - 5 . Together with the rela- 
tions N h * = N 1 (pN) h - 1 and N > N h * , the last term of 
Eq. I|A8() is approximated as 



1 + N/N h . ~ N/N^ =s e^(pN) 



(A9) 



Because 0.5 < L - h* - 0.5 < 1 from Eq. l|A"Tjl. we may 
further approximate it as 



ffr,\r\ L - h *-°- 5 ~ c t 



e?{pN) 



e 1 {l+pN) 



L-h*-Q.b 



C, (A10) 



where C is a small number of order (pN) L ~ h Sub- 
stituting Eq. IjAlOfl into Eq. (| AS|> and neglecting the sec- 
ond term C which is much smaller than the first term, 
we get 



k ci ~ e>7V(i + p N) 



£-2.5 



(All) 



This entrainment threshold agrees well with numerical 
data without any fitting parameter. Dividing Eq. I|A11I) 
by Eq. (JT^J, we obtain e 7 (l + pN)~ - 5 = c, which is 
about 0.67, 0.54 and 0.39 respectively for pN — 6, 10 and 
20. These agree well to c = 0.67, 0.60 and 0.45 obtained 
by numerical fitting (see Fig. 01. 



APPENDIX B: ENTRAINMENT THRESHOLD 
FOR BIDIRECTIONAL NETWORKS 

So far, we have considered only directed networks. 
However, similar exponential dependences are found also 
for bidirectional networks. In this section, we provide 
estimation of the entrainment threshold for bidirectional 
ER networks. Such networks are generated as follows. 
Independently for any i and j < i, we set Aij = 1 with 
probability p and Aij = otherwise, and put Aji = Aij. 
Self-connections are forbidden, so that An = 0. 

The entrainment threshold is obtained in a similar way 
as in Sec. IIVI We apply a global tree approximation in 
forward connections (which implies that the pattern of 
backward connections also take the same structure). Ev- 
ery oscillator then receives only one forward connection. 
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By assuming pN 3> 1, every oscillator (except those in 
the last shell) receives approximately pN backward con- 
nections [^3. Therefore, the only difference from the ER 
directed networks is that backward connections come not 
from the last shell but from the next shell. Thus, we get 
the following equations: 



/zr(0i -t) + kT(9 1 



1. 



— r(9 h -e h ^)+KT(9 h -e h+1 ) 

pN 



^r(fc - fe_o = 1. 



(Bl) 

1 for 2 < h< L, 
(B2) 

(B3) 



Using the approximations available for large pN, we ob- 
tain 



r(0 ft 



pN{(pN) 



L-h+l 



1} 



k{ P N - 1) 



for h > 2. (B4) 



The entrainmcnt thresholds is thus 



pN{{pN) 



L-l 



1} 



pN - 1 



(pN) 



L-l 



(B5) 



Stability analysis can also be performed in a same manner 
as in Sec. IIV CI We then get the same relaxation time 
as Eq. (|26|l . It is thus found that the dependences of 
the entrainment threshold and the relaxation time are 
also exponential in the bidirectional networks, but their 
functional forms arc slightly different from those for the 
directed networks. 
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